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It has been shown in our previous publication (jBlawzdziewicz et all 120031 ) that high- 
viscosity drops in two dimensional linear creeping flows with a nonzero vorticity compo- 
nent may have two stable stationary states. One state corresponds to a nearly spherical, 
compact drop stabilized primarily by rotation, and the other to an elongated drop stabi- 
lized primarily by capillary forces. Here we explore consequences of the drop bistability 
for the dynamics of highly viscous drops. Using both boundary-integral simulations and 
small-deformation theory we show that a quasi-static change of the flow vorticity gives 
rise to a hysteretic response of the drop shape, with rapid changes between the compact 
and elongated solutions at critical values of the vorticity. In flows with sinusoidal tempo- 
ral variation of the vorticity we find chaotic drop dynamics in response to the periodic 
forcing. A cascade of period-doubling bifurcations is found to be directly responsible for 
the transition to chaos. In random flows we obtain a bimodal drop-length distribution. 
Some analogies with the dynamics of macromolecules and vesicles are pointed out. 



1. Introduction 

Investigations of dynamical properties of fluid-fluid dispersions (e.g., emulsions dBorwankar fc Case 
1997[ lMasonl[l999T l and polymer blends (|Tucker III fc Moldenaerll2002t IWindhab etal. 
20051 )) require detailed understanding of the behavior of viscous drops in creeping flows. 



Such understanding is also crucial in developm e nt of new drop-ba s ed microfluidic system s 



( Whitesides fc Stroockl [20(51 iTan et all l200i ISong et all [2006i iGrigoriev efaH hoOfjf ) 



Therefore, dr op dynamics at small Reynolds numbers have been extensively studied ex- 
perimentally dTorza. Cox fc Masonll972;lBentrev fc Lealll986t[Bigio. Marks fc Calabresei 
1998tlGuido. Minale fc Maffettonefl2000t Cristini. Guido. Alfani. Blawzdziewicz fc Loewenberg 
200361: iGuido. Grosso fc Maffettond 120041). computationall y jRallison fc Acrivod 11978 



Kennedy et al.lll994ICristini et aZ.lll998l;IZinchenko et allu^lcristini et aLll200llj20036l: 

Renardvll2006j) and theoretically ( Barthes-Biesel fc Acrivosll973HRallisonll980tlBlawzdziewicz et al, 
2002L 120031 : IVlahovska et oi.ll2005l h 



These investigations revealed complex nonlinear drop dynamics resulting from the 
coupling between the drop shape and fluid flow. Examples of nonlinear phenomena that 
occur under creeping-flo w conditions include formation of self-similar nec k regions during 
a drop breakup process (jBlawzdziewicz et aZ.lll997t iLister fc Stondll998h . universal slow 
evolution of drops near the critical flow strength above which there are no stationary drop 



2 Young, Blawzdziewicz et al. 

shapes (jBlawzdziewicz et aLlll99 iBlawzdziewicz et aZ.ll2002T l , and existence 

of two branches of stable stationar y shapes of highly viscous drops in two-dimensional 

Stokes flows with nonzero vorticity (jBlawzdziewicz et al. 2003 ). 

As revealed by the analysis presented by Blawzdziewicz et al. ( 20031 ). there exists a 
flow-parameter range where a high-viscosity drop can either adopt a nearly spherical 
shape stabilized primarily by the rotational flow component or an elongated shape sta- 
bilized primarily by capillary forces. Abrupt changes of the drop shape from one state 
to the other can be used in manipulation of emulsion micro-structure and for control- 
ling the behavior of highly viscous drops in microfluidic devices. Due to discontinuous 
changes of emulsion microstructure, the bistable drop behavior may also significantly 
affect emulsion rheology. Moreover, the mechanism of the bi stability is of fundamental 
interest because of close analogie s to the dynamics of ve sicles (Misbah | 2006t|Mader et all 
Iviahovska fe Gracial [20071 ) and macromolecules (|Blawzd*z!ewical2006l ) in external 



flows. 

Due to the fundamental significance and because of potential applications (such as 
those mentioned above) it is important to explore the dynamics of highly vi scous drops in 



linear flows with nonzero vorticity. However, the investigations presented bv lBlawzdziewicz et al 



(2003) were limited to stationary drop shapes and stationary external flows. In the present 
study we focus on drop behavior in time-dependent flows. We elucidate the physical mech- 
anism that give rise to the bistable drop behavior and examine the consequence of these 
mechanism for drop response to time variation of the fluid vorticity. 

Th e system dynamics is i nvestigated via direct boundary-integral simulations ( PozrikidTsl 



1992t ICristini et alj|2001i iB lawzdzi cwiczl I2006D and by using a small-de formation ap 



proach lVlahovskal (|2003f k Iviahovska. Blawzdziewicz fc Loewenberd (|2005l ). In particular 
we show that the small-deformation equations with only several essential terms retained 
reproduce complex dynamical features of drop evolution that are associated with drop 
bistability. 

To emphasize important aspects of the drop dynamics we consider three flow varia- 
tion protocols. In the first protocol, the vorticity is slowly increased and then decreased. 
We find that such quasistatic vorticity ramping gives access to both the elongated and 
compact, nearly spherical stationary drop shapes. The drop exhibits a hysteretic behav- 
ior, with transitions between the compact shape (rotationally stabilized) and elongated 
shape (stabilized by capillary forces) occurring at different values of the vorticity when 
it is slowly ramped up or down. 

In the second protocol, the vorticity undergoes finite-frequency periodic oscillations. 
As expected, at low frequencies the drop behavior is quasistatic, with a hysteresis loop 
analogous to the one observed for linear ramping. At high frequencies the oscillations 
average out, and the drop undergoes small oscillations around the stationary shape cor- 
responding to the average flow. However, at intermediate frequencies wc find a much 
more complex behavior. In particular we show that there exists a frequency and am- 
plitude domain where the drop response to the periodic forcing is chaotic. Since in the 
creeping flow regime fluid motion is governed by the linear Stokes equations, the nonlin- 
ear chaotic behavior of the drop stems entirely from the coupling of the drop shape to 
the fluid velocity. An analysis of drop motion in the chaotic domain indicates that the 
transition to chaos is associated with the existence of two stationary states observed in 
steady flow. 

In our third flow- variation protocol, the vorticity of the imposed flow undergoes random 
changes. We observe, that the resulting statistical distribution of the drop length is 
bimodal in a certain regime of flow parameters, with two peaks around the drop length 
corresponding to the short and long stationary solutions. Hence, we find that also in 
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Figure 1. Decomposition of a linear incident flow into pure strain and rigid-body rotation. 

this problem the existence of two stationary states underpins drop behavior in a time- 
dependent flow. 

In §[2] the system considered in our paper is defined. The quasistatic hysteretic drop 
behavior is analyzed in §[3l our results for chaotic drop dynamics are presented in §[4j 
and drop motion in linear flows with randomly varying vorticity is discussed in §[5] Our 
findings are summarized in §[6l 



2. Viscous drops in two-dimensional linear flows 

We consider a viscous drop suspended in an incompressible Newtonian fluid of a con- 
stant viscosity \i. The viscosity of the drop fluid is fi = A/i, and the interfacial tension 
between the two phases is a. The drop is surfactant free, and no Marangoni stresses are 
present. There are also no buoyancy forces. We focus here on nonlinear effects that stem 
entirely from the coupling of the fluid flow to the drop shape (but not from the fluid 
inertia). Therefore, the creeping-flow conditions are assumed. 

In the creeping-flow regime the fluid motion in the regions inside (/ij = fj) and outside 
(//, = (J,) the drop is governed by the Stokes equations 

^V 2 u = Vp, (2.1) 
V-u = 0. (2.2) 

The fluid velocity u is continuous at the drop interface fl. Due to the absence of the 
Marangoni stresses the tangential viscous traction is also continuous. The jump in the 
normal viscous traction across fl is equal to the capillary pressure 

[n • r • n] = 2Ka, (2.3) 

where t is the viscous stress tensor, n is the outward normal unit vector, and k is the 
local curvature. 

The drop is subject to two-dimensional linear incident flow 

u (r)= 7 (E s +/3rj)-r, (2.4) 

where 7 is the strain rate, (3 is the dimensionless vorticity parameter, r is the position, and 
E s and f2 are the symmetric and antisymmetric parts of the normalized velocity-gradient 
tensor. In an appropriately adopted coordinate system we have 




E s = - 1 , n = - -1 , (2.5) 
z \ / z \ / 

without a loss of generality. Accordingly, (3 = corresponds to a purely straining flow with 
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Figure 2. Hysteretic evolution of viscous drop in two-dimensional straining flow with slowly 
varying vorticity. Viscosity ratio A = 200 and capillary number Ca = 0.20. (a) Normalized 
drop length I and (b) drop angle <j) versus vorticity parameter /3. Arrows indicate the direction 
of increasing time, and inset shows dependence of f3 on time. (Results from boundary-integral 
simulations.) 



the extcnsional axis x = y and the comprcssional axis x = —y, and /3 = 1 corresponds to 
shear flow in the x direction with the velocity gradient in the y direction. The tensor fi 
in equation (|2.4[) describes the rigid-body rotation in the anti-clockwise direction with 
the angular velocity 

w = \py. (2.6) 

The decomposition of the flow field (|2.4p into the straining and rotational components 
(|2.5p is sketched in figure [T] 

The dynamics of our system is characterized by three dimensionless parameters: There 
is the viscosity ratio A that characterizes dissipative forces in the drop- and continuous- 
phase fluids. The capillary number 

Ca=^i (2.7) 
(j 

(where a is the radius of an undeformed drop) characterizes the ratio between the de- 
forming viscous forces produced by the imposed flow (|2.4j) and the capillary forces that 
resist drop deformation and drive the drop towards the equilibrium spherical shape. Fi- 
nally the vorticity parameter (3 characterizes the magnitude of the rotational component 
of the external flow relative to the extensional component. 

In this paper we focus on the parameter regime where the drop deformation may be 
significant, which requires that Ca = 0(1). (However, the flow is not strong enough to 
cause drop breakup.) We also assume that the drop-phase fluid is much more viscous 
than the continuous-phase fluid, 

A»l. (2.8) 

Drop deformation process is hindered at large drop-phase viscosities, while drop rotation 
is only weakly affected by the viscous stresses inside the drops. Therefore the relative 
effect of the drop rotation is amplified in the regime (|2.8p . and the rotational component 
of the external flow produces nontrivial qualitative effects. 
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3. Hysteretic drop behavior 

3.1. Capillary and rotational stabilizing mechanisms 

To illustrate the effect of the rotational component of the flow (|2.4[) on the dynamics of a 
highly viscous drop we consider a system where the parameter (3 is first slowly increased 
and then slowly decreased. We adopt here a linear ramping protocol where the vorticity 
is slowly ramped up from (3 = to (3 = 0.4 and then ramped down back to zero. Before 
the ramping occurs the flow is maintained at = (for 5 % of the total ramping time) 
to allow the drop to relax to the stationary shape in purely straining flow. This vorticity 
variation protocol is represented in the inset of figure [2j a). 

The evolution of the drop shape in this time-dependent flow is depicted in figure [2l 
Figurc[2][a) represents the drop length I, and figurc[2j6) shows the drop angle 4> (measured 
anticlockwise from the axis a;) ; both quantities arc plotted versus the vorticity parameter 
j3. In our example, the drop viscosity is A = 200. The capillary number Ca = 0.2 is below 
the critical value for drop breakup (Ca = 0.22 in 2d straining flow) but it is sufficiently 
large to allow for a significant flow-induced drop deformation. The total ramping time 
is T = 2000t 7 , so that the drop response to the flow variation is nearly quasistatic. The 
calculations were performed using the adaptive boundary-integral procedure developed 



bv lCristini et all (|200ll . 119981 ). 

The results shown in figure indicate that the drop response to the vorticity variation 
is hysteretic. At (3 = drop is elongated, and it is aligned with the cxtensional axis of the 
straining component of the flow (</> = 7r/4). With increasing (3, the drop orientation slowly 
changes towards the symmetry axis x (0 = 0), and the drop length slowly decreases. At 
a critical value of the vorticity parameter, fii ~ 0.29, a discontinuous change occurs: the 
drop length and the angle suddenly decrease. Afterwards the drop is almost spherical and 
nearly aligned with the axis x. When the direction of the vorticity change is reversed, 
the drop initially retraces its trajectory. However, the drop does not jump back to the 
elongated shape until the vorticity reaches the lower critical value (3\ « 0.22 < /?2- 

The bistable drop behavior and the associated hysteretic shape evolution stem from 
the existence of two mechanisms that can stabilize a viscous drop in linear flows with 
rotation (|2.4[) . Namely, the drop can be stabilized by the capillary stresses (which drive 
the drop towards the equilibrium spherical shape) or by the vorticity flow component 
(which rotates the drop out of the cxtensional axis of the straining component of the 
flow). 

In a purely strain ing flow the drop assumes the interfacial-tension-stabilized elongated 
shape (|Tavlorlll934 ). The form of the drop results from the balance between drop defor- 
mation by the flow and drop relaxation due to the capillary forces. The deformation and 
relaxation occur on the respective time scales 

ty = \r\ (3.i) 

t a = X/jaa^ 1 , (3-2) 

both of which are proportional the viscosity ratio A ^> 1. The drop deformation D = 
(I — 2a) I a is determined by the time scale ratio 

D ~ tv/L, = Ca, (3.3) 

and therefore it is independent of the viscosity ratio in the limit A — > oo. In purely 
staining flows, the drop is oriented along the cxtensional axis x = y. 

For small values of (3, the vorticity flow component produces an 0{(3) perturbation of 
the drop orientation. The corresponding decrease of the drop length is 0((3 2 ). However, 
a further drop rotation is arrested because the straining component of the flow produces 



6 



Young, Blawzdziewicz et al. 



(a) (b) 




Figure 3. Schematic representation of the physical mechanism leading to bistable drop behavior 
in two-dimensional linear flows with non-zero vorticity. (a) Elongated drop is stabilized by 
capillary forces and destabilized by flow rotation; (6) compact drop is stabilized by flow rotation 
and destabilized by extensional flow component. 

hydrodynamic stresses that pull the elongated drop back towards the straining axis (as 
illustrated in figure EJa) . 

Since the 0(A _1 ) internal circulation inside an elongated high- viscosity drop is weak, 
the drop in its stationary state behaves analogously to a rigid body whose equilibrium 
orientation results from the balance of the torques produced by the straining and rota- 
tional components of the external flow. The transition to the compact drop shape occurs 
when the vorticity flow component becomes too strong to be balanced. Under such con- 
ditions, a rigid body would undergo a transition to a periodic motion with continuous 
rotation in the clockwise direction. Similarly, a drop also starts to continuously rotate 
when /3 achieves the upper critical value /?2- However, during the rotation the drop length 
decreases because the drop becomes misaligned with the extensional axis of the flow. As 
a result, the drop relaxes to a nearly spherical shape. 

In this new, compact shape the fluid inside the drop circulates with the angular ve- 
locity u>d that is nearly equal to the angular velocity of the external flow (|2.6|) . Within 
each period of rotation the drop undergoes a small deformation produced by the strain- 
ing component of the external flow (as schematically illustrated in figure [3^). However, 
the deformation does not accumulate because it is constantly convected away by the 
rotational component of the flow. 

Since the rotation occurs on the time scale 

t m t = (3.4) 

and the drop deforms on the much longer timescale p.ip . we find that the drop defor- 
mation 

d ~ w*7 = (w 1 ( 3 - 5 ) 

is small for A ^> 1, consistent with the results shown in figure ^a). 

Relation (|3.5p indicates that the deformation of the rotationally stabilized drop in- 
creases with the decreasing parameter /3. When (3 falls below the lower critical value /3i, 
the hydrodynamic torque associated with the straining component of the flow acting on 
a slightly elongated drop becomes strong enough to reorient the drop along the straining 
axis and arrest further drop rotation. Deformation thus starts to accumulate, the drop 
is stretched, and a transition to the intcrfacial-tcnsion stabilized elongated state takes 
place. 
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As shown in figure [2j a drop in the compact, rotationally stabilized stationary state is 
nearly aligned with the symmetry axis of the applied flow x. This behavior stems from 
the flow-reflection symmetry of Stokes equations and the fact that the drop is stabilized 
by rotation rather than the capillary forces. In the absence of the capillary forces (or 
in the limit of infinitely strong flow) the symmetry of Stokes equations implies that the 
stationary drop shape is invariant with respect to flow reflection. Hence, the shape is 
also invariant with respect to the corresponding transformation (x, y, z) — > (—x,y,z) 
of the spatial coordinates (and this symmetry corresponds to drop alignment in the x 
direction). A perturbation due to the presence of the capillary stresses produces only a 
small asymmetry because the effect of capillary forces is insignificant for a nearly spherical 
drop. 
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3.2. Parameter dependence of drop response 

The quasistatic response of the drop length to the variation of the vorticity is illustrated 
in figure 0] for different values of the capillary number Ca and viscosity ratio A. The 
results indicate that the size of the hysteresis loop is the largest for large values of A and 
Ca. When the capillary number is decreased, the upper critical vorticity parameter 02 
(corresponding to the transition from the elongated to the compact drop) decreases but 
the lower critical parameter (3\ remains nearly unaffected. In contrast, the viscosity ratio 
A affects primarily the position of the lower critical parameter (3\ (corresponding to the 
transition from the compact to the elongated drop). 

This behavior is consistent with the scaling relations (|3.3p and (|3.5p and the mecha- 
nism of drop bistability explained in § 13.11 According to these relations the deformation 
of an elongated drop scales with Ca and is independent of A, while the deformation of a 
compact drop is independent of Ca and scales with A -1 . We recall that the critical states 
corresponding to the transitions between the elongated and compact drop shapes corre- 
spond to the points where the maximal torque r 7 exerted by the extensional component 
of the flow on the drop just before the transition marginally balances the torque r ro t 
exerted by the vorticity flow component. Since r 7 ~ D whereas r ro t is approximately in- 
dependent of D, we conclude that /3i varies with A -1 and (3% with the Ca. This conclusion 
is consistent with our simulation results. 

The plots shown in figure [?] indicate that for a given viscosity ratio, there exists a 
critical value of the capillary number Ca* below which the drop response to the changes 
of the flow vorticity does not exhibit a hysteretic loop. The bifurcation point occurs at 
the critical value of the vorticity parameter (3* that corresponds to the position of the 
infinitesimal hysteresis loo p for Ca slightly above the c ritical value Ca*. It can be shown 
from the results derived by iBlawzdziewicz et al. ( 2003 ) that 



15 / 5\ 1/2 _ 3 /15 x 1/2 



for A 3> 1. The plots shown in figure[3]are consistent with the above expressions. We note 
that for A = 50 (see figure 0]c) the hysteretic drop behavior is not observed. The reason is 
that in this viscosity range Ca* exceeds the critical value of the capillary parameter for 
drop breakup. Hence, the elongated stationary drop shape does not exist for Ca > Ca*. 

3.3. Small-deformation analysis 

3.3.1. Evolution equations 

Crucial features of the evolution of a highly viscous drop in two-dimensional flows 
with nonzero vorticity are captured by small-deformation equations. In our approach 



(jBIawzdziewicz et aLll2003|) , the position rs of the drop interface is expanded in spherical 



harmonics 



r s /a = a' + V2^[// m Rc(y im ) + f^Tm(Y lm )}, (3.7) 

where I > and I > m > denote the spherical-harmonic order, f' lm and are the 
expansion coefficients, and the parameter a' is given by the drop-volume constraint. 
Since for m = all spherical harmonics arc real, wc set / ( q = in expansion (|3.7p . 
Moreover, flow-induced drop deformation preserves the symmetry of the incident flow 
(|2.4p . Therefore only even values of I and m need to be included in the analysis. 

Evolution equation for the expansion coefficients // and f" m are obtained by insert- 
ing the series (|3.7p into the boundary- value problem (|2.ip - (|2.5p . performing a boundary- 
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perturbation analysis, and reexpanding resulting products of spherical harmonics using 
appropriate Clebsch-Gordan coupling coefficients. The detailed analysis and explicit ex- 
pressions for the evolution equations at d ifferent truncation levels are presented elsewhere 
( Vlahovska|[200llVlahovska et a*Jl2005() . 



For simplicity, our small-deformation calculations are performed with the expansion 
(|3.7[) truncated at the lowest spherical- harmonic order 1 = 2, which leaves us with three 
independent drop-shape components: f 22 , f 22 , an d f 2 o- Noting that 

Re(Y~ 22 ) ~ cos 20, Im(y 22 ) ~ sin 20, (3.8) 

we find that the shape parameters / 22 and / 22 correspond to the drop deformation along 
the symmetry axis x and the straining axis x = y, respectively. The parameter / 20 
described an axisymmctric deformation along the axis z. 

The evolution equations for the shape parameters / 22 ,/ 22 , and / 20 , truncated at the 
second-order in the drop deformation, can be represented in the following form 

.4 = \-\d n + d 12 /^o)/22 ~ A^Ca" 1 [Drf^ - D 2 (f£ - fg - f£)] , (3.9a) 
f\ 2 = -2wf% 2 + A" 1 [^1/22/22 - Ca-^Dx + 2D 2 f^)f 22 ] , (3.96) 
f>> 2 = 2uf! 22 + A" 1 [(d 31 + d 32 / 20 + d^f! 2 2 + d 34 /22 + ^3 5 / 22 ) - Ca- 1 ^ + 2D 2 f^)f^ 2 ] , 

(3.9c) 

where the dot denotes the time derivative (normalized by 7 -1 ). The terms involving 
the coefficients dy correspond to drop deformation by the external flow, and the terms 
involving Dk describe the capillary relaxation. All these terms are 0(A _1 ) in the large- 
viscosity-ratio re gime. Explicit expressions for the coefficients dij and Dk are given in 



( Vlahovskall2003l ) : here we only note that these coefficients are functions of the viscosity 
ratio A and have finite limits for A — > 00. 

The two remaining terms on the right-hand side of equations (|3 . 9 £>[) and (|3.9cj) (the 
terms proportional to oj) are viscosity independent. These terms represent the rigid-body 
rotation of the drop, with the angular velocity 

a, = -§j8+| Cl /£ a , (3.10) 

where c\ = (15/271") 1 / 2 . Consistent with our qualitative physical picture described in ij l3.ll 
(and illustrated in figure [3|), the rotational velocity (|3.10|) involves two terms. The first 
term corresponds to the rotation of the drop by the vorticity component of the flow 
(|2.4p . The second term, which is proportional to the shape parameter f 22 that described 
deformation in the x direction, corresponds to the rotation of a deformed drop by the 
straining component of the external flow towards the straining axis x = y. 

3.3.2. Reduced description 

It has been shown by iBlawzdziewicz et ~d\ (120031) that for A ^> 1 the drop behavior 



near the bifurcation point p.6p can be described by simplified asymptotic equations 

f' 22 = -2w# a - A^Ca-^Di/M, (3-lla) 
/ 2 ' 2 = 2^/ 22 -A- 1 Ca- 1 J D 1 / 2 ' 2 + A- 1 J 3 i. (3.116) 

where D\ = 20/19 and J31 = (57r/6) 1//2 are the high-viscosity limits of D\ and d^\. The 
asymptotic result (|3.11[) is obtained form (|3.9p on assumption that near the bifurcation 
point there is a balance between drop deformation and rotation (which corresponds to 
A -1 ~ w/22) an d the balance between capillary relaxation and rotation (which yields 
(ACa) _1 / 22 ~ / 2 ' 2 and (ACa) _1 / 22 ~ / 22 ). Moreover, the two contributions to the angular 



10 



Young, Blawzdziewicz et al. 



velocity (|3.10|) are of the same order but do not cancel (i.e., u> ~ j3 ~ f^)- Equations 
(|3.9p are rescaled accordingly, and only the leading order terms are retained. 

Equations (|3 . 1 1 [) have all necessary ingredients that are needed to describe the hys- 
teretic drop behavior. We have terms representing drop rotation by the straining and 
vorticity components of the external flow (i.e., the terms proportional to to), drop relax- 
ation due to the presence of the capillary forces (the terms proportional to Ca -1 ), and 
stretching of the drop along the straining axis of the straining component of the external 
flow (the last term in equation (|3 . 1 1 £>[) ) . Neglecting any of these terms would qualitatively 
alter the solution structure, which no longer would manifest the key features of the drop 
evolution in the parameter range considered herein. 

3.3.3. Asymptotic solution 

In the regime A -1 < 1 and /?< 1 the stationary solutions of equations (|3.11[) can be 
obtained by a singular-perturbation analysis. To the leading order in the small parame- 
ters, we find that the elongated drop is described by the relations 

w~0, (3.12a) 

A _1 Ca _1 Di/^ ~ A _1 d 3 i. (3.126) 

The first of the above relations correspond to the fact that an elongated drop does not 
rotate (there is only a weak fluid circulation inside it, as predicted by our qualitative 
analysis). The second relation describes the balance between drop deformation by the 
external flow and relaxation due to the capillary forces. Recalling the definition (|3.10p of 
the angular velocity u> we find that relations (|3.12|) yield 

/aa - ci/9, (3.13a) 

/& ~ D^dn Ca. (3.136) 

By inserting the above relations back into (|3.11[) one can verify that they constitute a 
consistent leading-order asymptotic stationary solution. 

According to equation (|3.136[) drop elongation along the straining axis x = y (i.e., 
= 7r/4) scales with the capillary number, which is consistent with the scaling result 
(|3~5|) . The drop angle 

= i arctaa(/^/^ a ) (3.14) 
only slightly deviates from <j> = it/ 4 because f 2 2 <^ f^2 (assuming that Ca ^> (3 and 
/?«1). 

The leading-order stationary solution corresponding to the compact drop is obtained 
from the following relations 

- 2w/& - 0, (3.15a) 

2^/22 — -A^efeij (3.156) 
which are obtained by dropping from the evolution equations (|3.11| the 0(X~ 1 ) capillary- 
relaxation terms. Taking into account the definition (|3.10| of w we thus obtain 

f'22 ~ (3.16a) 

and 

(the solution with the plus sign in front of the square root is unstable). 



Hysteretic and chaotic dynamics of viscous drops 



11 



1 .5 
1 .4 

1 .3 

1 .2 



(a) 
Ca - 



-0.20 

0.18 
0.16 

0.14* 

0.10 



.0 
0.0 





0.2 




0.4 




Figure 5. Quasistatic variation of drop length I with vorticity parameter (3 for A = 200 and dif- 
ferent values of capillary number (as labeled), (a) Solution of small-deformation equations (|3.9p ; 
(b) comparison of small-deformation results (dashed lines) with boundary integral simulations 
(solid lines). 



Since the shape parameter f 22 vanishes according to equation (|3.16q|1 . the drop is 
oriented in the x direction. For A -1 <C p wc find 



which is consistent with our scaling estimate (|3.5[) . For p < pi, where 

/3 1 = 2(c 1 J 3 i) 1/2 A- 1 / 2 , 



(3.17) 



(3.18) 



the solution (|3.166[) does not exist; the drop thus undergoes a transition to the long 
solution. The critical vorticity parameter decreases with increasing A, in agreement with 
our numerical results presented in figure |4] 

The solution (|3. 12j) - (|3. 18[) of the simplified small-deformation equations (13.111) is per- 

turbative. We note, however, the exact stationary solution can also be found (jBlawzdziewicz et al 
20031) . Our analytical solution quantitatively agrees with the results of numerical simu- 
lations, provided that the drop deformation is not too large. 

3.3.4. Numerical results 

Predictions of the small-deformation equations (|3.9[) for drop behavior in a two-dimensional 
linear flow with slowly varying vorticity are depicted in figure [5] for a system with the 
viscosity ratio A = 200. Figure^ a) shows the dependence of the drop length 



l/a 



(3.19) 



on the vorticity parameter p for the same set of capillary numbers as those represented in 
figure [U Figure [5f 6) compares the small deformation results directly with the results of 
the boundary-integral simulations. The small-deformation calculations were performed 
using the second-order equations (|3.9p because for an elongated drop they are more 
accurate than the simplified equations (|3.12[) . 

The results shown in figures |4] and [5] indicate that for small and moderate capillary 
numbers the small-deformation theory yields accurate quantitative predictions. At high 
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Figure 6. Evolution of drop length (left) and angle (right) in two dimensional linear flow with 
harmonic variation of vorticity <|4.ip . for different values of period T normalized by drop-de- 
formation time (as labeled). Mean vorticity (5 — 0.25, vorticity amplitude 5/3 = 0.13, viscosity 
ratio A = 275 and capillary number Ca = 0.2. Panels (c) depict chaotic dynamics. (Results from 
mall-deformation theory.) 



values of Ca drop behavior is also captured quantitatively, except for the upper portion 
of the hysteresis loop (i.e. when the drop is in the elongated state). For all values of 
the capillary number the lower and upper critical vorticity parameters f3\ and fii are 
obtained within the numerical error of the boundary- integral simulations. Our additional 
calculations (not shown) indicate that a similar accuracy is obtained for A = 50 and 
A = 100. 



4. Chaotic drop dynamics in a sinusoidal straining flow 

Dynamical systems with multiple equilibrium states often exhibit novel dynamics when 
driven by simple forcing ( Guckenheimer fc Holmea 1983a). Thus, despite the laminar 
nature of the Stokes flow, we expect to find interesting nonlinear dynamics of a viscous 
drop in a time-varying linear flow with rotation. To explore this dynamics we will now 
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Figure 7. Period doubling in the dynamics of viscous drop in two-dimensional linear flow with 
harmonic variation of vorticity. Vorticity parameter (3 (solid lines), drop deformation D (dashed), 
and drop angle <f> (dash-dotted) are shown versus time t normalized by the oscillation period T. 
Viscosity ratio A = 275, capillary number Ca = 0.2, period T/t 7 = 1.14, mean vorticity f3 = 0.21, 
and vorticity oscillation amplitude (a) 8(3 = 0.6 and (6) 0.8. (Results from small-deformation 
theory.) 



investigate the drop response to harmonic variation of the vorticity 

0{t) = (3 + 5(3 cos(2irt/T), 



(4.1) 



where (3 is the average vorticity value, 6(3 is the oscillation amplitude, and T is the 
oscillation period. 

We have performed a series of small deformation calculations (§ 14. 1[) and boundary- 
integral simulations (§ I4.2|) for different values of the flow parameters (3, 8(3 : and T. If the 
oscillation period T is much shorter than the drop deformation and oscillation times (|3.1[) 
and (|3.2p . we find that the drop undergoes small oscillations about a stationary shape 
corresponding to the mean value of (3 (which is an expected behavior). In the opposite 
regime T 3> tj,t ai the quasistatic drop behavior described in §[3] is recovered. In what 
follows we focus on the most interesting parameter domain T ~ tj,t a and 5/3 ~ (3\ — (32, 
in which an interaction of different timescales as well as an interplay between the short 
and elongated drop shapes is anticipated. 



4.1. Small-deformation results 

Figure[H]illustrates the dependence of the drop evolution in linear flow with the oscillatory 
vorticity (|4.1[) on the oscillation period T. The viscosity ratio in this example is A = 275, 
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Figure 8. Bifurcation diagram showing period doubling cascades and transition to chaos for 
viscous drop in linear flow with harmonic vorticity variation with mean (3 = 0.21 and period 
T/t 1 = 1.14. Viscosity ratio and capillary number are the same as in figures |B] and [7] (Results 
from small-deformation theory.) 

and the capillary number is Ca = 0.2. The mean value of the vorticity (3 = 0.21 is close 
to the lower critical value (3\ = 0.18, and the oscillation amplitude is 8(3 = 0.13. 

Figure^ a) represent our results for the shortest oscillation period of the flow vorticity 
T/tj = 0.36. The drop oscillates about the compact stationary shape in this case. Both 
the drop lengths and the drop angle vary periodically, with the period T<j equal to the 
period T of the external forcing. The drop length decreases when the drop is in the 
compressional quadrant — tt/2 < cf> < and increases for < <f> < tt/2. Upon an increase of 
the period of the external forcing the amplitude of the angular drop oscillations increases. 
When the oscillation amplitude reaches tt/2 the drop starts to tumble, as illustrated in 
figure b) . 

A further increase of the period of the external forcing results in a qualitative change 
of drop response. We find that the drop still undergoes a tumbling motion; however, the 
evolution is not periodic but it becomes chaotic, as shown in figure [HU c) . The chaotic 
motion continues up to T/t 7 = 1.14, and then the drop reverts to periodic motion. For 
T/tj (figure ®l) the drop oscillates about the elongated stationary shape, and in the 
regime T/t 7 3> 1 (figure [6)2) the system approaches the quasistatic behavior discussed in 

§eEI 

The transition to the chaotic drop motion occurs through a cascade of period doubling 
events, as illustrated in figures [7| and [H] Figure [7J depicts the drop evolution at three 
values of the amplitude of vorticity oscillations 8(3 (the remaining system parameters are 

f The only significant deviation from the quasistatic evolution for T/t 7 3> 1 occurs right after 
the drop jumps from the long to the compact shape when (3 increases above the upper critical 
value /?2- Namely, before the drop settles down to the compact stationary shape it undergoes a 
tumbling motion with the amplitude decaying on the timescale (|3.2p . 
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Figure 9. Same as figure |6j except that the results are from boundary-integral simulations 
and for slightly different values of oscillation period. Same (3 and 5(3 as in figure [5] 



the same as those that yield the chaotic motion depicted in figure [S]c). The results shown 
in figure [7£ a) indicate that for a sufficiently small oscillation amplitude {5(3 = 0.06 in 
our example) the drop evolves with the period 7d = T (i.e., the period equal to that 
of the external forcing). At a larger amplitude S(3 — 0.08 the drop oscillation period 
is T d = 2T and for 8(3 = 0.083 we find T d = AT. The period-doubling scenario of the 
transition to chaos in our system is further supported by the bifurcation diagram shown 
in figure [3 where the drop length I at times t = nT (n = 1,2,...) is plotted versus the 
flow-oscillation amplitude 8(3. 

An analysis of the results shown in figure [7] indicates that the period doubling occurs as 
a result of a resonance between the drop tumbling motion and the vorticity oscillations. 
Namely, if the drop is relatively long and approximately aligned with the straining axis of 
the external flow when the vorticity parameter (3{t) reaches a minimum, the drop rotation 
may be significantly slowed down or even arrested (as in the long-drop stationary state 
discussed in §[3j) . Such a temporary arrest of drop rotation corresponds to the shoulders 
in the plots of the angular evolution depicted in figure [71 On the other hand, if the 
drop angle exceeds 4> = tt/2 when the vorticity goes through a minimum, the arrest of 
drop rotation does not occur. As seen in figure [7j b) this interplay of drop tumbling with 
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Figure 10. Bifurcation diagram for viscous drop in linear flow with harmonic vorticity variation 
with mean f3 — 0.22, T/t 7 = 1.82, A = 275, and Ca = 0.2. (Results from boundary-integral 
simulations.) 



oscillations of the external forcing produces the period-doubling bifurcation that leads 
to alternating accelerated and retarded drop-rotation cycles. 

In our numerical calculations depicted in figures [B][5] we have used the full set of the 
small-deformation equations (|3.9p but we find that the simplified asymptotic equations 
(|3.1ip yield similar results. In particular these equations correctly reproduce the cascade 
of the period doubling bifurcations and the chaotic-evolution domain. However, a further 
simplification of the evolution equations is not possible: if any of the terms in equations 
(|3.11| is removed the solutions qualitatively change and the chaotic domain disappears. 
We also find that chaotic drop dynamics occurs only for highly viscous drops with A > 
200. 

4.2. Boundary-integral results 

Figurc[5]shows examples of our simulation results for the periodic and chaotic drop evolu- 
tion. The flow parameters are similar to those used in our small-deformation calculations 
described in ij l4.ll As with the small-deformation calculations, for short periods of the 
external forcing T the drop oscillates around the compact stationary shape, for moderate 
periods the system undergoes a transition to chaotic evolution, and for long periods the 
drop motion approaches the quasistatic behavior. Consistent with the small-deformation 
results, the chaotic drop dynamics revealed by the boundary-integral simulations is due 
to a cascade of period doubling bifurcations. A bifurcation diagram illustrating this be- 
havior is presented in figure [TU] 

The domain of chaotic motion found in the direct boundary-integral simulations some- 
what differs from the corresponding domain obtained from the small-deformation the- 
ory. Also, the magnitude of the chaotic fluctuations in the drop length is larger in the 
boundary-integral runs. We have tested the convergence of the boundary-integral simu- 
lations, and we believe that the differences between the drop behavior obtained by the 
two different methods stem from the approximations involved in the small-deformation 
theory. We note, however, that the evolution in the chaotic and period-doubling regimes 
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Figure 11. Length probability distribution for viscous drop with A = 200 and Ca = 0.2, in 
linear flow with stochastic vorticity. Mean value of the vorticity parameter (j3) = 0.25, variance 
S/3 = 0.13, and correlation time r corr /t 7 = 0.24. Inset shows vorticity probability distribution. 
(Results from boundary-integral simulations.) 

is very sensitive to the to small changes of the system parameters (hence, also sensitive 
to the approximations involved in our calculations). 




5. Drop Statistics in a linear flow with Stochastic Vorticity 



In some systems (e.g. emulsion flows through a packed bed of fibers (jMosler fc Shaqfeh 



19971) or turbulent emulsion flows with drops that are much smaller then the Kolmogorov 



scale ( Cristini et al. 2003 of )) a viscous drop undergoes deformation in a random external 
creeping flow. To gain some understanding of the role that drop bistability may play in 
such systems we consider the drop behavior in a flow with stochastic vorticity. 

We assume that the time variation of the vorticity parameter (3 is described by a sta- 
tionary Markovian Gaussian process (i.e, the Ornstcin- Uhlenbcck process) wit h the mean 



), variance 8/3, and correlation time r corr . A standard numerical scheme ([Fox et 



19881 ) for generating such a time-correlated Gaussian process is used to model the time 



variation of (3 for a given drop trajectory. 

An example of drop behavior in a stochastic flow with a Gaussian variation of the 
vorticity is presented in figure [TTJ The mean value of the vorticity is ((3) = 0.25 and the 
variance is 5(3 = 0.13. The correlation time of the vorticity distribution t coxt /t a = 0.24 
is several times shorter than the drop relaxation time. The figure depicts the probability 
distribution of the drop length, obtained using the boundary-integral simulations. The 
capillary number is Ca = 0.2, and the drop viscosity is A = 200. 

The results indicate that the drop-length distribution is bimodal for the above param- 
eter values. This behavior is expected since the vorticity undergoes random variation 
in the domain that includes both the lower and upper critical values j3i, and 02 of the 
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Figure 12. (a) Length and (b) angle probability distributions for viscous drop with A = 200 and 
Ca = 0.2, in linear flow with stochastic vorticity. Mean and variance of the vorticity parameter 
are the smae as in figure [TT1 the correlation times are r corr /t 7 = 0.025, 0.1, 0.2, 0.3, 0.4 for lines 
marked 1-5, respectively. (Results from small-deformation theory.) 



vorticity parameter (3. Since drop response to a flow with slow variation of vorticity is 
hysteretic, a drop in the random flow tends to stay in the neighborhood of the compact 
and the elongated stationary states. 

We note that the peak of the length probability distribution at I sa 1.25 is shifted 
towards the shorter drop lengths compared to the length of a drop in the elongated 
stationary shape. This is because Ca = 0.2 is close to the critical capillary number for 
drop breakup. Th us, due to the slow time s cale in the drop dynamics near the critical 
capillary number ijBlawzdziewicz et al\ 120021 ). the drop does not have sufficient time to 
fully extend before the vorticity significantly changes. 

Figure [T^] shows the probability density distribution for the drop length I and drop 
angle <$ for different values of the flow correlation time r corr . Other system parameters are 
the same as in figure 1111 The calculations were performed using the small-deformation 
equations (|3.9|) . The results indicate that at short flow correlation times the drop-length 
probability distribution is peaked around small values corresponding to the short-drop 
stationary solution, and the has a moderate- hight peak at 4> ~ it /A. As the flow correlation 
time T corr increases, the length probability distribution becomes bimodal: one of its peaks 
corresponds to compact and the other to elongated drops. A corresponding change occurs 
in the angle distribution, i.e., its peak becomes more pronounced, and shifts towards the 
straining axis <fi = it/ A. 

The shift of the typical drop length and orientation from the compact to elongated 
state when the flow correlation time is increased resembles the analogous shift for a 
system with harmonic vorticity oscillations (sec figures [6^ and figures [6]ff) . This behavior 
is further illustrated in figure Q1J] which shows the average values and the variance of I 
and cb versus the correlation time r corr . 



6. Conclusions 

We have presented results of numerical and theoretical investigations of the dynamics 
of highly viscous drops in two-dim ensional linear creeping fl ows with time-dependent 
vorticity. In our earlier publication (|BIawzdziewicz et all 120031 ) we predicted that in sta- 
tionary flows such drops exhibit bistable behavior: there is a range of system parameters 
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Figure 13. (a) Mean and (b) variance of drop length (dashed) and angle (broken line) versus 
correlation time T corr normalized by drop deformation time, for A = 200, Ca = 0.2, (0) = 0.25, 
and 5/3 = 0.13. 



where the drop may assume either an elongated shape approximately aligned with the 
straining axis of the external flow or a nearly spherical shape approximately aligned in 
the flow direction. Here we analyze the consequences of this behavior for the system 
dynamics. We also elucidate the physical mechanism that leads to drop bistability. 

A direct consequence of the existence of two stationary states is hysteretic drop re- 
sponse to a flow with slowly varying vorticity. We have explained that the rapid transition 
from an elongated non-rotating drop shape to the nearly spherical compact shape occurs 
when the vorticity becomes strong enough to overcome the effect of the straining flow 
component that aligns the drop with the straining axis. This transition is thus analogous 
to the behavior of an elongat ed rigid parti cle which starts to tumble when the vorticity 
grows above a critical value ( Jefferv 1922( ). A viscous drop also begins to tumble at a 
critical vorticity magnitude fa. However, when the drop becomes misaligned with the 
extensional axis it relaxes under the action of capillary forces towards a nearly spherical 
rotationally stabilized stationary shape. 

If, in turn, the vorticity is slowly decreased, the drop returns to the elongated shape 
only after the vorticity magnitude reaches a lower critical value fa < fa- This hysteretic 
drop response occurs at high drop viscosities because the rotational stabilizing mecha- 
nism is more efficient in the high- viscosity regime. A highly viscous drop deforms less 
within each drop revolution, so the compact shape remains stable even for small vorticity 
magnitudes. 

The existence of two stationary states affects drop dynamics not only in the quasistatic 
regime but also at finite frequencies of the external forcing. At small amplitudes of har- 
monic vorticity oscillations the drop simply oscillates (with the same frequency as the 
external forcing) about one of the stationary states. However, if the vorticity-variation 
range includes both critical values fa and fa the dynamics of the system is much richer. 
We find that with an increasing magnitude of the vorticity oscillations the system under- 
goes a cascade of period-doubling bifurcations resulting in chaotic drop dynamics. The 
period doubling stems from the resonance between the periodicity of the external forcing 
and the tumbling motion of the drop when it jumps from a (partially) elongated shape 
towards the compact rotationally stabilized state. 

Chaos in our system emerges despite linearity of Stokes equations - the system dy- 
namics is nonlinear because of the coupling of the flow to the evolving fluid interface. A 
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detailed analysis of small-deformation equations describing drop dynamics reveals that 
in addition to chaos associated with the period-doubling mechanism there also exists 
in our system a different kind of chaotic evolution that results from manifold tangling 
( Guckcn heimer fc Holmes! [l98 3fr). Our analysis of different types of chaos will be pre- 



sented in a separate publication. 

To our knowledge, chaotic drop dynamics in Stokcs-flow regi me has never been observed 



before. We note, however, that in a recent independent study (jKas-Danouche et aZ.ll2007f> 
has reported chaotic dynamics in co-annular Stokes flow with insoluble surfactant ad- 
sorbed on the fluid interface. Chaos in their system also appears as a result of period 
doubling. 

Understanding of drop bistability and the associated dynamical phenomena is relevant 
for many practical problems. For example, interpretation of rheological response of emul- 
sions of highly viscous drops to time-varying flows requires insight into drop dynamics. 
Our results may also be useful in design new methods for manipulating emulsion mi- 
crostructure in material processing and controlling drop behavior in microfluidic flows. 
Drop bistability could, e.g., be used to construct microfluidic switches, and chaotic drop 
dynamics may be relevant for microfluidic mixing. 

The stabilizing and destabilizing mechanisms described in our paper apply not only 
to viscous drops but also to other deformable particles. Therefore, results of our study 
have a broader significance. 

In particular, our analysis suggests that macromolecules with high degree of internal 
dissipation may undergo a transition between a nearly spherical and moder ately elon- 
gated states (in addition to the standard coil-stretch transition predicted by Ide Gennesl 
1974) . In fact, the dynamics of macromolecules can be modeled using equations analogous 



to (|3.11[) . supplemented with terms representing random thermal forces. Such a simpli 
fied description correctly captures the most important features of power s pectra of DNA 



mole cules evolving in linear flows with nonzero rotational component (jBlawzdziewicz 
20061) . 



There are also close analogies between drop and vesicle motion. The main difference 
between these two systems is that vesicles satisfy a constant-area constraint whereas 
the drop area can vary. This constraint gives rise t o periodic vesicle motion (such as 



tank treading and tumbling) even in stationary flows (|Misbahll2006l ; IVlahovska fc Gracia 



20071) . It would be interesting to determine if a coupling of vesicle oscillations to a har- 
monic variation of the external flow can lead to chaotic dynamics. 

It would also be of significant interest to experimentally explore the bistable and chaotic 
drop dynamics (as well as related phenomena for other deformable particles). In such 
experiments a four-roll mill could be used to produce a linear flow with a controlled 
magnitude of vorticity. The experiments could al so be performed using recently develo ped 



microfluidic analogues of a four-roll mill device ([Hudson et aZ.ll2007t iLee et alluOOlv 
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